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Abstract 

A discrete solid-on-solid model of epitaxial growth is introduced which, in 
a simple manner, takes into account the effect of an Ehrlich-Schwoebel barrier 
at step edges as well as the local relaxation of incoming particles. Furthermore 
a fast step edge diffusion is included in 2 + 1 dimensions. The model exhibits 
the formation of pyramid-like structures with a well-defined constant inclina- 
tion angle. Two regimes can be distinguished clearly: in an initial phase (I) 
a definite slope is selected while the number of pyramids remains unchanged. 
Then a coarsening process (II) is observed which decreases the number of is- 
lands according to a power law in time. Simulations support self-affine scaling 
of the growing surface in both regimes. The roughness exponent is a = 1 in 
all cases. For growth in 1 -|- 1 dimensions we obtain dynamic exponents z = 2 
(I) and z = 3 (II). Simulations for d = 2 seem to be consistent with z = 2 (1) 
and z = 2.3 (II) respectively. 



Among the different morphologies of growing surfaces in molecular beam epitaxy, 
the frequently observed formation and coarsening of pyramid-like structures with a 
well-defined slope is of particular interest [0, |^, ||, H, If. This kind of behavior 
has been seen in experiments of molecular beam epitaxy in such diverse systems as 
Fe(OOl) [I, g, Cu(OOl) g, GaAs(OOl) g |10|, and HgTe(OOl) The selection 

of a constant inclination angle can be associated with the compensation of slope- 
dependent uphill and downhill currents of deposited particles. This has been studied 
by numerical integration of continuum equations |^, ^ . 

In the framework of our discrete model, particles are deposited on a (i-dimensional 
simple cubic lattice with L'^ sites. The prescription obeys a solid-on-solid restriction, 
i.e. neither overhangs nor holes will occur and the growing structure is completely 
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characterized by integer height variables assigned to each of the sites. The time t 
associated with the growth process is defined as t = N/L'^ where is the total 
number of deposited particles. We will consider a relaxation mechanism for the 
incoming particles as well as an additional diffusion process with high energy barriers 
at step edges. Only after the particle has reached its final position, a new one is 
deposited on the surface. 

To begin with, the model is explained in terms of a one-dimensional substrate 
{d = 1) which is initially fiat {hj{t = 0) = for j = 1, 2, . . . , L) and satisfies periodic 
boundary conditions. A new particle is placed with equal probability anywhere on 
the surface, say at site i. If any of the neighbor sites j (with \j — i\ < Rinc) offer a 
lower height hj < hi, the freshly deposited particle is moved to the lowest of these 
sites (a random one in case of a tie). This models the effect of particles reaching the 
substrate with a residual momentum and kinetic energy |12]. The parameter Rinc is 
called the incorporation radius of the model. 

Without additional diffusion this prescription would coincide with the well-known 
Random Deposition with Surface Relaxation (RDSR) 0. RDSR belongs to the 
Edwards-Wilkinson universality class the scaling behavior of which is known ana- 
lytically 1^, 13, 0]. In particular, the surface becomes flat in 2 + 1 dimensions, i.e. 
its width increases only logarithmically with time and system size. 

However, if the incoming particle is not yet bound to in-plane neighbors after 
relaxation, it is considered to diffuse along the surface. Now, an energy barrier at step 
edges can effectively surpress downward moves since the above mentioned residual 
momentum has been lost. In the following, this Ehrlich-Schwoebel barrier is taken 
to be inflnite. Consequently, the diffusion is restricted to connected sites of equal 
height hi. The main result of this Letter is that the interplay of this diffusion and 
the incorporation mechanism yields a stable slope which is controlled by Rinc only. 

The second parameter of the model is the diffusion length which marks the 
maximum distance that the particle can explore before it comes to rest. A simplistic 
prescription is used to replace the actual diffusion in 1 + 1 dimensions: if a potential 
in-plane neighbor is available within a distance the particle is placed there, other- 
wise it sticks to the deposition site. This mimics the nucleation of diffusing particles 
with a mean free path of order l^j |T5[, a quantity which depends on the ratio of the 
diffusion constant and the incoming flux (see e.g. [0). 

The diffusion process alone (without relaxation) does not alter the frequencies at 
which values hj occur in comparison with simple Random Deposition (RD). Thus, 
formally, ( /i ) = t and w{t) = {{h^)-{h ff/"^ oc with {3 = 1/2 where (...) de- 
notes the average over all lattice sites. The width w of the surface grows indeflnitely 
and no saturation due to lateral correlations is observed. However, the spatial dis- 
tribution of /i-values is changed drastically. The inflnite Ehrlich-Schwoebel barrier 
results in the formation of separate islands of linear size oc on the flat substrate. 
Without local relaxation no coarsening is possible, the number of islands and their 
sizes remain constant after a short transient. However, their height and thus the 
slopes grow indeflnitely resulting in the RD power law. The proflle of a single island 
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Figure 1: Growth in 1 + 1 dimensions: evolution of a single island for various system 
sizes L = Id = lOO(o), 200(n), 400(A), and 800(+). Results were averaged over 
50 independent runs, error bars would be smaller than the symbols, a) shows the 
average local slope m vs. scaled time for three different values of the incorporation 
radius (top to bottom: Rinc = 1, 2, 3). In b) the evolution of the surface width w is 
displayed for the same system sizes and Rmc = 1- 



is given by the inverse of a cumulative binomial distribution (i.e. an inverse error 
function in the continuum limit) as was observed independently by Krug |T^. 

The net effect of the infinite Ehrlich-Schwoebel barrier is an uphill current of 
deposited particles which can be compensated for by the relaxation process described 
above. The resulting slope selection is most clearly studied for Id = L when only 
one island forms. Figure 1(a) shows for this case the evolution of the average local 
height difference m = ( | hj^i — hj\) for three different values of Rmc- (The profiles 
produced by our model are strictly monotonic along the slopes of the pyramids, 
otherwise m would systematically over-estimate their inclination.) The curves for 
different system sizes L = Id collapse when time is scaled with indicating that 
after a characteristic time oc /|, a constant slope is reached. 

Note that if a particle is deposited on a terrace of length ra, it jumps downward 
with probablity Rinc/n whereas it is added to the upper terrace with probability 
1 — Rinc/n. Therefore, the selected slope for which uphill and downhill currents 
cancel is rrisat = {'^Rinc)~^ in excellent agreement with the simulations. A more 
detailed presentation of this argument, also in terms of a continuous description of 
the step dynamics, will be given in a forthcoming publication [Q. 

While the single island forms, the width w of the corresponding surface increases 
until it reaches its saturation value Wsat = T^sat L / a/48. Thus, the roughness ex- 
ponent is a = 1 as expected for a constant slope. Simulations suggest a growth 
exponent (3 = 1/2 which corresponds to a dynamic exponent z = 2 in agreement 
with the time scale found for the slope selection process. Apparently, the relaxation 
process is basically irrelevant for the small time behavior and w increases as in Ran- 
dom Deposition. Figure 1(b) displays the evolution of w for various system sizes, 
rescaled according to the assumed behavior. 
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Figure 2: The dynamics of the model for /^i = 20 ^ L in 1 + 1 dimensions: a) surface 
width vs. time, rescaled according to the assumption z = 3 and a = 1 for system 
sizes L = 200(n), 283(x), 400(A), 566(0), and 800(+). Averages were performed 
over 50 runs with standard errors smaller than the symbol size. In b) the evolution 
of w{0), lO-m(x), and Uisi (A) is shown for a system with L = 5000 and Id = 20 on 
average over 40 simulations. After the slope is selected, surface width and number 
of islands evolve according to w oc t^^^ and riisi oc 



Next we investigate the dynamics of many islands. For 1 <K lo <<C L the number 
of initially formed pyramids is of the order of L/l^. The small time behavior is 
dominated by the slope selection process which is completed after the characteristic 
time tm oc l"^. In the second phase fluctuations in the deposition of particles induce 
a coarsening process in which smaller islands are merged into larger ones [0. The 
number of islands riisi decreases according to a power law: Uisi oc This 
exponent coincides with theoretical predictions for noise driven coarsening in one- 
dimensional systems p, |T9| and can be related to domain growth in magnetic systems 
2l| , Accordingly, the surface width increases like w oc t^^'^. After a time of 



order the system reaches its saturation state: a single remaining pyramid with 
the above given constant inclination. 

Fig 2(a) shows results of simulations with = 20. Under a rescaling according to 
z = 3 and a = 1 the curves for different system sizes collapse for times greater than 
the L-independent value tm- The above scaling law is very well confirmed by Monte 
Carlo results provided the number of islands is still fairly large. Fig 2(b) shows the 
results for L = 5000 and In = 20 with riisi decreasing from about 150 to 10. A more 
detailed analysis indicates that, close to saturation, the pairwise merging of the last 
few islands occurs at characteristic times and yields a cascade-like increase of w. 

It is important to note that the condition Id ^ 1 has to be satisfied in order 
to produce the observed behavior. Otherwise, nucleation is likely to occur even on 
small terraces, the above mentioned monotonicity of the profile is lost, and riisi can 
fluctuate strongly in time. The detailed behavior of systems with a small diffusion 
length Id = 0{1) will be the objective of forthcoming investigations. 

In the following we discuss a generalization of the model to 2+1 dimensions. We 
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assume Rmc = 1 and the diffusion process is simulated explicitely as a random walk 
under the constraint that the particle remains at the same height. The walk stops 
when the particle reaches an in-plane neighbor or has accomplished ud jumps. This 
corresponds to a typical diffusion length Id ^ ^/nD provided the set of accessible 
sites has a compact shape (unlike, for example, a fractal cluster). 

Hereupon the particle performs an additional diffusion (random walk at constant 
height) along the step edge until it is bound to (at least) two in-plane neighbors or 
until moves have been accomplished. We would like to emphasize that without 
step edge diffusion the properties of the model are changed drastically. For instance, 
we find a roughness exponent a ~ 0.8 for = whereas the expected slope 
selection with a = 1 is confirmed for n£)2 ^ 1. The crossover is in itself interesting 
but here we concentrate on the model with l£)2 oc ^nm large enough (of order V) 
to ensure that the particle reaches the nearest accessible kink site. In this regime 
the bases of the growing pyramids are virtually of square shape, cf. Fig. 3, whereas 
smaller values of Im yield rounded cross-sections. 

The overall behavior is qualitatively the same as in 1+1 dimensions. First, a 
constant slope builds up which is observed most clearly in the formation of a single 
pyramid (y^n^ V). Figure 4(a) shows that here the surface width obeys a scaling 
relation w = ■ g{t/L'^) where g{x 2> 1) = const. Again, slope selection occurs 
(a = 1). The systematic deviation from the prediction nisat = 1/2 is of order 1/L in 
finite systems, cf. inset in Fig. 4 (b). This correction is taken into account in Fig. 4 
by scaling w with {rrisatL) rather than with L. 

The observed dynamic exponent is in good agreement with z = 2. li g{x) (x 
for X -C 1, this would imply a/ z = j3 = 1/2 as in 1 + 1 dimensions. Yet, the simple 
power law seems to hold only for very short times. This can be seen clearly in the 
inset of Fig. 4(a) where the effective exponent obtained from two consecutive data 
points is plotted. As already observed in 1+1 dimensions, the surface width increases 
initially like in RD with (3 = 1/2. The relevance of the incorporation mechanism is 
determined by the number of step edges which is independent of the pyramid size 
in 1+1 dimensions. However, for d = 2 it increases with the island radius already in 
the submonolayer regime. Thus, the deviation from RD-behavior becomes noticeable 
much earlier. 

Preliminary results concerning the coarsening of many islands {Id ^ L) are 
shown in Fig. 4(b). The data seem to be consistent with a dynamic exponent z ?a 
2.3 {j3 ~ 0.43) in this regime. Simulations with small and intermediate values of 
nD2 indicate that (3 is significantly smaller when the step edge diffusion is slower. 
However, a more precise determination of this value is necessary and should be 
based on the simulation of larger systems. At the current stage a definite conclusion 
concerning the scaling behavior cannot be drawn. The above mentioned analogy 
to the dynamics of domains in magnetic systems would imply 1/ (3 = z = ?> but is 
unlikely to carry over to epitaxial growth in 2+1 dimensions ||^. Several authors 
obtain values (3 < 0.25 from the numerical studies of continuum equations (see 
and refs. therein). A recent argument by Tang et ai. suggests that z = 2 + d for 
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Figure 3: Snapshots (density plots) of the evolution of the surface morpholgy in 
d = 2 at times t = 4, 16, 64. The system size is L = 100, the model parameters are 
Rinc = 1,Id = 15, and Id2 = 200. 



noise-driven coarsening. 

The relevance of our simple model for molecular beam epitaxy is readily rec- 
ognized if one identifies the parameter with the typical island distance in the 
submonolayer regime. Several authors have calculated this length for the simple 
case in which two adatoms already form a stable nucleus (||16| and refs. therein). 
For = 2 one obtains h (x {D / Ff {{D / Ff in d = I H) where D is the diffusion 
constant and F the flux of arriving particles. A helpful simplification of our model is 
due to the infinite Ehrlich-Schwoebel barrier. However, even with a finite (but non- 
vanishing) barrier the qualitative scenario will be the same. Extensions of the BCF 
theory |2^ which include an incorporation mechanism also yield slope selection 
There, however, the tops of pyramids are rounded. Furthermore the terrace 



width depends on the actual energy barriers and the temperature of the system 

In summary, we have presented a simple model which reproduces essential fea- 
tures of epitaxial growth. Local relaxation leads together with surface diffusion in 
the presence of a high Ehrlich-Schwoebel barrier to the selection of a stable slope 
of growing islands. We identify two different scaling regimes (both with roughness 
exponent a = 1): the initial slope selection (I) and the coarsening of pyramids (II). 
For growth on a one-dimensional substrate we find z = 2 in regime (I) and 2; = 3 in 
(II). Simulations in d = 2 suggest the same dynamic exponent for (I) and a value of 
z 2.3 in regime (II). Apart from a more accurate measurement of these exponents 
in larger systems, forthcoming investigations will address the crossover behavior for 
Id = C'(l), the precise influence of the step edge diffusion length Id2, the effect of a 
finite Schwoebel barrier, and a continuous description of the step dynamics |18 



We would like to thank S. Kohler and M. Schroder for useful discussions. This work 
has been supported by the Deutsche Forschungsgemeinschaft through SFB 410. 
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Figure 4: Scaling in 2+1 dimensions for system sizes L — 140(+), lOO(o), 70(n), 
and 50(0) with 1^2 = 2L. Results were averaged over 30 independent runs with 
error bars smaller than the symbols. Fig. a) shows the surface width vs. time for 
one island {^/nD ~ i^), rescaled according to z = 2.0 and a = 1. The inset displays 
the effective exponent calculated from two consecutive data points with L = 50. In 
b) the surface width is plotted vs. time for the model with coarsening {no = 200), 
rescaled according to z — 2.3 and a — 1. The inset shows the saturation slope vs. 
1/L. 
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